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The transport of heavy, neutral or charged, point-like particles by incompressible, resistive mag- 
netohydrodynamic (MHD) turbulence is investigated by means of high-resolution numerical simu- 
lations. The spatial distribution of such impurities is observed to display strong deviations from 
homogeneity, both at dissipative and inertial range scales. Neutral particles tend to cluster in the 
vicinity of coherent vortex sheets due to their viscous drag with the flow, leading to the simultaneous 
presence of very concentrated and almost empty regions. The signature of clustering is different for 
charged particles. These exhibit in addition to the drag the Lorentz-force. The regions of spatial 
inhomogeneities change due to attractive and repulsive vortex sheets. While small charges increase 
clustering, larger charges have a reverse effect. 

PACS numbers: 52.30.-q, 52.65.-y, 52.30.Cv 



I. INTRODUCTION 

Many natural and teehnologieal settings involve the 
transport of particles with mass densities much higher 
than the main constituents of the underlying turbulent 
plasma. For instance in astrophysics, dense dust grains 
are suspended in the interstellar medium^ or in molec- 
ular clouds.'^ Despite usually constituting only a small 
fraction of matter in a given astrophysical system, dust is 
of central importance for an understanding of many pro- 
cesses, such as the heating and cooling of interstellar and 
intergalactic plasmas.— It is also of significance for the 
detection of interstellar magnetic and velocity fields and 
their turbulent properties i^i^ Dust grains also play impor- 
tant roles in turbulent protoplanetary disks for the for- 
mation of planetSj?- for the interaction of comets with the 
solar wind)^ or the modification of Kolmogorov spectra of 
weakly turbulent shear Alfven waves in the ionospherci^ 

In industrial applications, fusion plasmas usually con- 
tain a significant population of heavy charged parti- 
cles mostly coming from the introduction of dust in the 
plasma by erosion of walls. ^^d^d^d'^ An understanding of 
dust dynamics in this context is not only of interest in 
order to quantify the effect of impurities but also to pos- 
sibly remove the dust particles by exploiting the notion 
that scattering through waves can enhance the drag force 
on dust particles by orders of magnitude. "'^^ 

An important aspect of the passive transport of dense 
particles is the fact that their inertia can rarely be ne- 
glected and is in general responsible for the apparition 
of strong inhomogeneities in their spatial distribution. 
Indeed, dense particles experience a centrifugal accelera- 
tion in the rotating regions of the flow and tend to con- 
centrate in those places where strain is dominant. In 
the purely hydrodynamical turbulent motion of neutral 



fluids, this well-studied phenomenon is usually referred 
to as preferential concentration.^''^ Its main manifestation 
is the creation of large particle voids in the most ex- 
cited parts of the turbulent flow, which are separated by 
high-concentration regions where the particles follow the 
calmer motions of the fluid. The process of ejection from 
eddies is universal and is also expected in the case of con- 
ducting fluids. It is well known, however, that in MHD 
turbulence, the large-strain regions are characterized by 
the presence of current /vorticity sheets where both the 
magnetic field and the fluid velocity experience strong 
fluctuations (see, e.g., Ref . [l^ . Consequently particles 




FIG. 1: (color online) Density of heavy particles in MHD 
turbulence: gray (blue online) : low density, black (red online) : 
high density. 

tend to cluster in the vicinity of these sheets, as already 
observed in anisotropic MHD simulations at low mag- 
netic Reynolds numbers. In statistically isotropic sit- 
uations, such an observation is confirmed as illustrated 
in Fig. [T] where sheet-like structures can clearly be ob- 
served in the instantaneous three-dimensional distribu- 
tion of neutral particles obtained from a direct numerical 
simulation. 
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These correlations between the most active regions of 
the flow and the particle distribution contrast much with 
what is observed in neutral flows where, conversely, the 
particles flee coherent structures. This implies major dif- 
ferences in the properties of dense particles suspensions 
depending whether the flow is charged or not. For in- 
stance, we will see that inertial-range particle dynamics 
in MHD is ruled by concentration processes rather than 
ejection processes, leading to larger fluctuations in their 
spatial mass distribution. 

This paper is organized as follows. Section [Til gives a 
detailed description of the model we consider, of its do- 
main of validity, and of the numerical methods we use 
to integrate it. In Section Ull] we present a qualitative 
comparison of the effects of preferential concentration 
in charged and neutral turbulent flows, emphasizing in 
each case the different roles played by the flow coherent 
structures. Section llVI is dedicated to a more quantita- 
tive study of particle clustering, both at dissipative and 
inertial-range scales. In the former case, we show that 
particles cluster on a singular set displaying multi-fractal 
properties. In contrast, this scale- invariance does not 
hold in the inertial range where the mass distribution 
rather depends on a scale-dependent contraction rate of 
the particle dynamics. In Section |Vl we present results 
on the effect of the particle charge on their concentration 
properties. Finally, Section |VT] encompasses concluding 
remarks. 



II. MODEL AND METHODS 

We consider an electrically conducting fluid whose dy- 
namics is described by the MHD equations 

dtu + u ■ Wu ^ X B) X B - + f + t^V^M, (1) 
5tB V X (m X B) + 7 + r/dV^B, (2) 
V • M = 0, V • B = 0, (3) 

where u is the fluid velocity field and B is the mag- 
netic field written with the dimension of a velocity (i.e. 
it is the magnetic field expressed in the Gaussian sys- 
tem of units divided by ^/Airpf where pi is the fluid mass 
density), v and ry^ are the kinematic viscosity and the 
magnetic diffusivity, respectively. / and 7 are two forces 
such that the amplitudes of the Fourier modes associated 
to wavelengths of modulus less than 2 are kept constant. 
Such a choice provides an input of kinetic and magnetic 
energies, preventing them to condensate at large length 
scales. This means in particular that there is no mean 
magnetic field, (B) = 0, so that the turbulence is sta- 
tistically isotropic. The electric field E is related to the 
magnetic field by the Ohm's law 



E = 7;dV X B u X B. 



(4) 



We consider next a spherical particle embedded in the 
MHD fiow described above. We assume that its radius 
a is much smaller than the smallest active scale of the 
carrying fluid flow (in turbulence, the Kolmogorov dissi- 
pation scale 77 = (j^^/ek)^^^ where ek designates the mean 
kinetic energy dissipation rate) . We also assume that the 
mass density pp of the particle is much higher than the 
fluid mass density pf. These assumptions allows one to 
approximate the particle by a point and considering that 
the main effect exerted on it by the fluid is a Stokes vis- 
cous drag, which is proportional to the velocity difference 
between the particle and the fluid flow. When in addi- 
tion the particle is uniformly charged with density pc and 
evolves in an electro-magnetic field, it is subject to the 
Lorentz force, so that its trajectory X{t) is a solution to 
the Newton equation 



X = 



1 



E{X,t) + XxB{X,t) , (5) 



Here E has the dimension of a square velocity, meaning 
in the Gaussian system, that it has been multiplied by 

c/y/Wpi. 



where the dots stand for time derivatives and t = 
2/9pa^/(9/C)fi/) and d = pp c/ {pcy^4:Trpf) (with the charge 
density pc expressed in Gaussian units). Note that we 
neglect here the effect of gravity. 

The first term on the right-hand side of ([5]) represents 
the viscous drag of the particle with the flow. It in- 
volves the Stokes time r, which is the relaxation time 
of the particle velocity to the fluid velocity. This char- 
acteristic time scale gives a measure of particle inertia 
and is usually non-dimensionalized by the Kolmogorov 
time scale ~ (jz/e^)^/^ to define the Stokes number 
St^r/r,^ = (2/9)(pp/pf)(a/r;)2. When St<^ 1, the par- 
ticle responds very quickly to the flow fluctuations and 
almost follows the dynamics of tracers, i.e. X k, u{X,t). 
Conversely, when St ^ 1, a long time is required for 
the particle to react and consequently it moves quasi- 
ballistically with a low velocity. 

The second term on the right-hand side of (O involves 
the length-scale £, which measures the importance of the 
Lorentz force exerted by the ambient fiuid magnetic and 
electric fields onto the particle. By analogy to the Stokes 
number which parameterizes the relative motion of par- 
ticles at scales of the order of the Kolmogorov dissipa- 
tive scale ?7, a non-dimensional Lorentz number can be 
introduced as Lo = t Brms/^ (where -Brms is the root- 
mean square (r.m.s.) value of the magnetic field). Note 
that, conversely to St, this number involves a large-scale 
quantity, iJrms, because of the cross-product quadratic 
non-linearity appearing in the Lorentz force. This num- 
ber can be written as Lo = (A^/tt /9) B^-^s Pcci^ / (ci^y/Pf) 
and measures the relative effect of the Lorentz force with 
respect to the viscous drag. When Lo <C 1, the Lorentz 
force is negligible and inertia dominates. When Lo 3> 1, 
the former is dominating the particle dynamics. 

Note that, besides the hypotheses of small radius and 
large mass density that are necessary to write the viscous 
drag in the present form, it is implicitly supposed that 
the particle charge does not affect the local motion of the 
charged fluid. For this, it is required that the Coulomb 
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100 


0.073 


0.14 


5.6 ■ 10""* 


8.6 • 10"* 


5-10-^ 


2.45 ■ 10"^ 


2.2 • 10"^ 


0.94 


2.7 


17 


256^ 
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TABLE I: Parameters of the numerical simulations. Re\ = y^l5V L/v. Taylor-Reynolds number, iirms: r.m.s. velocity, Brms: 
r.m.s. magnetic field, tk: mean kinetic energy dissipation rate, em'- mean magnetic energy dissipation rate, v. kinematic 
viscosity, rjd- diffusivity, Sx: grid-spacing, 77 — (i^'^/tk)^''*: Kolmogorov dissipation length scale, = (z//ek)^^^: Kolmogorov 
time scale, L = (2/3i5)^/'^/(ek -I- em): integral scale, = L/uj-ms- large-eddy turnover time, A^'^: number of collocation points, 
Np-. number of particles of each species. 



force exerted by the particle onto the plasma charges that 
are at a distance of the order of the particle radius is much 
smaller than the Lorentz force exerted on these charges 
by the ambient electromagnetic field. Namely, one must 
have (4/3)7rpca ^ (1/c) Wrms-Srms \/47''/Of- This condition 
restricts the model validity to not-too-high values of the 
particle charge density, and thus of the Lorentz number 
Lo. However, some freedom is left for the choice of pa- 
rameters. Indeed, the charged particle is characterized by 
three parameters: its radius a, its mass density Pp and 
its charge density pc and one can show that it is always 
possible to find a combination of these parameters satis- 
fying the three constraints on the validity of this model 
with arbitrary values of the Stokes and Lorentz numbers. 

Finally, in the proposed model, we neglect the collec- 
tive effects of the particles onto the MHD flow, i.e. the 
particles are passive. This implies that the application 
fields are restricted to very diluted situations where the 
particle number density is so small that their hydrody- 
namical and electromagnetic effects on the conducting 
carrier fluid can be disregarded. This does not prevent 
from considering, for instance for statistical purposes, a 
large number of such passive particles in the flow. 

The numerical simulations are carried out using the 
reformulation of the momentum equation ^ in terms of 
the vorticity u) — Wxu and a pseudo-spectral solver. The 
underlying equations are treated in Fourier space, while 
convolutions arising from non-linear terms are computed 
in real space. A Fast-Fourier- Transformation (FFT) is 
used to switch between these two spaces. The time 
scheme is a Runge-Kutta scheme of third order. The 
inter-process communication uses the Message Passing 
Interface (MPI). The main parameters of the simulation 
are given in Tab. HI 

In addition to the fluid fields, we introduce three dif- 
ferent types of particles into the flow: fluid tracers cor- 
responding to the degenerate case St = Lo = 0, heavy 
uncharged particles with Lo = and 11 different val- 
ues of St between 0.05 and 20, and finally heavy charged 
particles with St = 1 and 8 different (positive) values 
of Lo between 0.015 and 1.5. One million particles of 
each species are uniformly seeded into the statistically 
stationary ffow. After a period of relaxation of 2.4 T^, 
where Tl desigantes the large-scale eddy turnover time, 
the statistical analysis of particle dynamics is started and 
averages are performed over a time span also equal to 
2.4 Tf,. In order to obtain from the grid values the veloc- 
ity and magnetic fields at the particle positions we use 



a tri-linear interpolation. This interpolation scheme par- 
allelizes efficiently and provides sufficient accuracy^ for 
the statistical properties under consideration. 

III. PREFERENTIAL CONCENTRATION IN 
CONDUCTING FLOW 

We focus in this Section on the effect of the Stokes drag 
on the qualitative properties of particle clustering. For 
this we consider uncharged particles (with Lo = 0), so 
that the Lorentz's force in the Newton equation ^ can 
be neglected. In an incompressible flow, the dynamics of 
tracer particles {St = 0), whose trajectories are solutions 
to X = u{X,t), is volume preserving. Consequently, 
an initially uniform distribution of tracers remains uni- 
form at any later time. Because of their drag with the 
fluid, heavy particles have a dissipative dynamics and 
their spatial distribution tends to develop strong spatial 
inhomogeneities, which are observed to span almost all 
scales of the fluid turbulent motion (see, e.g., Fig. [T]). 
The particles strongly deviate from the fluid motion in 
those places where the latter experiences spatial or tem- 
poral fluctuations. The clustering properties are thus 
essentially coupled to the statistical properties of veloc- 
ity gradients (and accelerations). At small scales (below 
the Kolmogorov scale 77 where the flow is differentiable) 
the spatial distribution of particles is settled by the one- 
point fluid gradient statistics. At inertial-range scales, it 
is the structure of the gradient field that determines the 
spatial organization of particles. Hence, depending on 
the observations scales, different processes are invoked to 
describe clustering. 

The small-scale clustering can be mostly explained us- 
ing the formalism developed for dissipative dynamical 
systems. As is well known trajectories of such systems 
tend to concentrate onto a fractal set called attractor 
(see, e.g., Ref. Hfll). For heavy particles whose trajec- 
tories solve a non-autonomous second-order differential 
equation, this set is embedded in the position-velocity 
phase space {X,X) and evolves dynamically. The in- 
stantaneous positions of particles are obtained by pro- 
jecting this set onto the position subspace. The fractal 
properties of the attractor can be preserved by this pro- 
jection, leading to non-trivial scaling laws for the parti- 
cle distribution at small scales. The fractal properties 
of the attractor are fully determined by the statistical 
properties of the linearized particle dynamics, which 
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only involves the fluid velocity gradient along particle 
paths. The conducting character of MHD turbulent flow 
does not alter such a mechanism. Although the prob- 
ability density function (PDF) of velocity gradients in 
MHD quantitatively differs from that observed in neutral 
fluid turbulence, the qualitative picture of dissipative- 
scale clustering is the same. For instance, in both cases, 
concentration processes are optimal for Stokes numbers 
of the order of unity. Such qualitative similarities are 
confirmed in next Section. 




FIG. 2: (color online) Isosurfaces of vorticity (online: green) 
and particles (online: red) in a slice in purely hydrodynamic 
turbulence (Left) and in magnetohydrodynamic turbulence 
(Right). In both cases, the domain has approximatively the 
size 256 ri x 256 and the slice width is ~ 5r; 

Strong qualitative differences between neutral and con- 
ducting carrier flows show up in the spatial distribution 
of particles at scales within the inertial range. In hy- 
drodynamical turbulence this distribution anti-correlates 
with the spatial organization of the rotating structures. 
As illustrated in Fig. [2] (Left) rotating structures where 
vorticity is strong tend to expel particles that then con- 
centrate in the calmer regions in between eddies. The 
mechanism of ejection from vortex filaments is sketched 
in Fig. [3) The flow structure of a conducting fluid is very 
different. In MHD turbulence violent gradients corre- 
spond to high strain regions and not to strong rotations. 
The coherent structures are thus vorticity/current sheets 
and the particles tend to agglomerate on them as illus- 
trated from Fig. [2] Right. As sketched in Fig. [3] a heavy 
particle tends to approach closer to a vortex sheet than 
fluid tracers, clearly leading to high concentration in its 
vicinity. Moreover, this phenomenon can be amplified 
when heavy particles cross the sheet. As vortex sheets 



almost correspond to discontinuities of the velocity field, 
heavy particles might hence enter a region with a veloc- 
ity very different from that of the fluid. This implies a 
phase during which the particle decelerates and which is 
responsible for the concentration of particles close to the 
sheets. Of course, the particle clusters and the sheets do 
not have a one-to-one correspondence. The fluid velocity 
close to the sheet is parallel to it and thus entrains the 
particles to the edges of the structure, a picture that ap- 
plies also to fluid tracers.— The fraction of particles that 
are close to a given current sheet depends on the dynam- 
ical history of the latter and in particular on its lifetime 
and on the fact whether or not it has experienced violent 
accelerations. 




FIG. 4: Vorticty second-order moment along particle paths 
normalized by the Eulerian moment, as a function of Stokes 
number. □: MHD turbulence; o: hydrodynamic turbulence. 



This different sampling of high-vorticity structures by 
heavy particles in hydrodynamic or MHD turbulence has 
important impact on the distribution of vorticity along 
their paths. This can clearly be observed by measuring 
vorticity fluctuations at the particle positions. Figure [3] 
shows the second-order moment of vorticity as a function 
of the particle Stokes number. While in hydrodynamic 
turbulence one observes for Stokes numbers of the order 
of unity a large drop-off of more than 30% with respect 
to the Eulerian vorticity, a gain of 15% can be measured 
in MHD turbulence. In both cases this trend displays a 
maximum for St of the order of unity and tends to dis- 
appear at large inertia. The effect of inertia on vorticity 
sampling is not visible in the core of the vorticity distri- 
bution. It is rather due to a change in the time fraction 
that the particles spend in the violent regions with large 
fluctuations of the fluid velocity gradient. 




IV. QUANTITATIVE STUDY OF CLUSTERING 



FIG. 3: Sketch of the motion of a heavy particle (bold line) 
in a vortex filament (left) and in a vorticity sheet (right); 
reference fluid tracer trajectories are drawn as thin lines. 



We now turn to more precise and systematic ways to 
characterize particle clustering. The particle spatial dis- 
tributions are of different natures if the fluid velocity fleld 
is differentiable or not. Dissipative and inertial-scale clus- 
tering are thus distinguished here. 
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A. Small-scale clustering 

At those scales where the fluid velocity field is differen- 
tiable, particle clustering stems from a competition be- 
tween the stretching imposed by the fluid velocity gradi- 
ent and the linear relaxation/dissipation due to the parti- 
cle viscous drag. This materializes by the convergence of 
particle trajectories to a dynamically evolving fractal at- 
tractor in the position-velocity phase space. The particle 
locations are then obtained by projecting such a singu- 
lar set onto the configuration space. Hence, depending 
on the fractal dimension, the projection is or is not also 
singular, leading to possibly fractal spatial distribution 
of particles. 

A (multi)fractal distribution is characterized by a di- 
mension spectrum Dp (see, e.g., Ref . [20|) . which is related 
to the scaling behavior of the various moments of the 
coarse-grained mass distribution. Consider the coarse- 
grained density pr defined as the mass of particles inside 
a fixed box of size r, divided by the volume of the box. 
The p-th order moment of the coarse-grained density be- 
haves as (pP) oc Ap~^'> i^p-'i) for r <^ rj, where d denotes 
the space dimension. Di is usually refered to as the infor- 
mation dimension and D2 as the correlation dimension. 
We focus here on second-order statistics, because they 
relate to the probability of finding two particles close to 
each other and, thus, to estimate the rate of local binary 
interactions (as e.g. collisions, chemcal reactions, etc.). 
Indeed, one can easily show that the probability P^{r) 
of finding two particles closer than a given distance r be- 
haves at small scales as P^{r) oc r^^ (see, e.g., Ref. [20|) . 
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FIG. 5: Correlation dimension for heavy, uncharged particles 
in MHD turbulence as a function of the Stokes number. □: 
average local slope of P2~{r), +: phase-space correlation di- 
mension D2 obtained from the fit P2~{r) — Ar^^+Br^. Inset: 
same for hydrodynamic turbulence (from Ref. l23h. 



Figure [5] represents the behavior of the correlation di- 
mension D2 as a function of the Stokes number for neu- 
tral (uncharged) particles. D2 is observed to decrease 
from d — 3 a,t St ~ 0, where the dynamics of parti- 
cles recover that of uniformly distributed tracers in an 
incompressible flow. It reaches a minimum for Si ~ 1 
(maximum of small-scale clutering), and increases again 
to reach 3 at large St where the dynamics of particles be- 



comes less and less influenced by the flow and approaches 
the ballistic motion of free particles. In order to estimate 
the correlation dimension, two different methods were 
used. The first one consists in averaging the local slope 
(i.e. the logarithmic derivative) of P^{r) over a scaling 
range of roughly two decades in r. The error is then esti- 
mated through the maximum and minimum deviations of 
the local value of the slope from its average. This method 
leads to over-estimate clustering for St ^ 1. Indeed, as 
was observed in synthetic flows^i the probability of flnd- 
ing in the physical space two particles at a distance less 
than r behaves as P^ (r) ~ Ar^'^ + Br"^, where D2 is the 
phase-space correlation dimension. When D2 < d, the 
physical-space correlation dimension is D2 — £'2- Con- 
versely, when D2 > d, the particle distribution is space 
fllling (to second order) and D2 = d. Assuming such a 
form for the two leading terms in P^ (r) gives much bet- 
ter estimates of the correlation dimension and allows one 
to observe a saturation of the physical-space correlation 
dimension to d = 3 for St > Stcr — 4. As seen from the 
inset of Fig. [5l the correlation dimension as a function of 
the Stokes number has a similar behavior as in purely hy- 
drodynamic turbulence. The main difference is a rescal- 
ing of the horizontal axis by a constant factor and comes 
from the deflnition of the Stokes number, which involves 
the turnover time associated to the Kolmogorov dissipa- 
tive scale. The measurement of the latter as the typical 
scale of the velocity gradient gives some uncertainty that 
might explain this discrepancy. 

The scaling behavior of the coarse-grained density dis- 
tribution, which is related to the dimension spectrum, 
implies a speciflc large-deviation form for the PDF of 
Pr- Indeed, one can easily see that if the moments 
behave as (p^) oc r'^"^^ '^p^'^', the local dimension 
hr = logpr/logr has to follow at small scales r the 
large deviation form p{hr) oc exp[— (logr) H(hr)], where 
H{hr) = maxp[(p— 1) {Dp — d)—phr] (see, e.g., Ref. I21I) . 
The rate function H is positive and convex and attains 
its minimum (equal to 0) for hr — Di. 



B. Inert ial-range clustering 

The large-deviation form for the distribution of the 
coarse-grained density that is discussed above gives a 
very precise meaning for the scaling behavior of the fluc- 
tuations of mass at small scales. Above the Kolmogorov 
scale, where the flow is itself almost self-similar and spans 
a large range of timescales, such an invariance is broken 
because particles tend to preferentially select those ed- 
dies of the flow whose associated turnover time is of the 
order of their response time. Figure [H] represents the var- 
ious PDFs of the coarse-grained density pr at fixing the 
particle Stokes number and varying the box size r inside 
the inertial range of the turbulent carrier flow. One ob- 
serves that the deviation from a Poisson distribution is 
maximal at small scales and that the particle distribu- 
tion approaches uniformity when r increases. It is clear 
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from this picture that the PDF associated to two differ- 
ent scales do not differ only through a scaling factor, so 
that there is no large deviation principle ruling the mass 
distribution in the inertial range. 
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FIG. 6: PDFs of the coarse-grained particle density over a 
scale r for neutral particles with Stokes numbers St — 1 and 
various values of r as labeled. 

The PDF of the coarse-grained density has noticeable 
characteristics that give insights on the mechanisms lead- 
ing to particle clustering. Both tails (at low and high 
values of the mass) are much wider than what would be 
expected from a uniform distribution. This implies that, 
on the one hand, it is pretty probable to observe a re- 
gion with a very large mass, but on the other hand it is 
also very probable to observe a region almost empty of 
particles. The large-density tails decrease exponentially 
or slower. This is in contrast to what was observed in 
hydrodynamic turbulence, where the decrease is clearly 
faster than exponential.^-'^^ This exhibits a difference in 
the clustering process and might be traced back to the 
distinct role of the coherent structures discussed in pre- 
vious section. The behavior of the small-density tails is 
far from being completely settled. As seen from Fig.[7l it 
seems to display a power-law decay, but the data seems 
spoiled by a plateau due to the limited number of par- 
ticles used in the simulation. For hydrodynamic turbu- 
lence where a larger number of particles were used, this 
effect was not so visible and clear evidence of a power-law 
behavior could be obtainedi^ 



As seen from Fig. [71 a uniform distribution is also re- 
covered if, in place of increasing r while keeping St fixed, 
one fixes r and decreases St. Strong deviations from uni- 
formity are indeed observed for 5f ~ I and the distri- 
bution approaches that of Poisson when St — *■ 0. This 
leads to asking the question whether or not the two lim- 
its of small inertia (r/ry) — > oo and r — s- are equivalent, 
or more precisely if the probability distribution p{pr) de- 
pends in this limit on a single parameter oc St/ (r/ry)". To 
answer this question, one first notices that in the limit of 
small inertia, the particle velocity is to leading order^*' 

^u{X,t)-T{dtu + u-Vu){X,t). (6) 

This approximation consists in assuming that the parti- 
cles behave as if they were advected by a synthetic com- 
pressible flow. The important time scale that measures 
the strength of particle inertia with response time r at a 
given scale r is given by the inverse of the rate 7(t, r) at 
which a particle blob of size r contracts when advected 
by this synthetic compressible flow, namely 



_ 1 



V • v d'^a; = T 



\x \ <r 



f V • (m • Vm) A^x. (7) 

J|a;|<r 



The idea is now to understand phenomenologically how 
the integral on the right-hand side behaves as a func- 
tion of r. One can first notice that using the diver- 
gence theorem, the volume integral can be rewritten as 
the surface integral on |a;| = r of w • Vu. We next as- 
sume that the blob overlaps a vortex/current sheet where 
the velocity almost experiences a jump. This implies a 
behavior of the surface integral ~ u-cmsiur) / t/) , where 
is the typical velocity difference over the Kolmogorov 
scale 77. In terms of the Stokes number St = t/t^, the 
blob contraction rate thus behave as 7 ~ Stuj-^s/i" and 
is thus proportional to the non-dimensional contraction 
rate P = {rj/r)St. Finally, this approaches predicts that 
the inertial-range distribution of particles in the limit 
St ^ Q or r /rj ^ 00 depends only on the non-dimenional 
contraction rate V = {ri/r)St <^ \. As seen from Fig. [H 
numerical results confirm this behavior. 



FIG. 

scale 
bers. 




7: PDFs of the coarse-grained particle density over a 
r — I877 for neutral particles with several Stokes num- 




FIG. 8: Probability distribution of the coarse-grained density 
for three different values of the contraction factor P = (rj/r) St 
(widest distribution: P ~ 0.07, intermediate P ~ 0.024, and 
narrowest: P ~ 0.006). Inset: deviations from uniformity of 
the variance of the coarse-grained density as a function of P. 
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V. EFFECT OF PARTICLE CHARGE 

Despite cases were impurities can be handled as neutral 
particles in a conducting fluid)^ there are several situa- 
tions where the particles have a charge density different 
from the carrier flowfii^ In this section we investigate the 
clustering properties of charged particles arising from a 
non-vanishing Lorentz number Lo. To shed light on the 
effect of a charge density different from the surrounding 
flow, we focus on ensembles of particles with identical 
Stokes number but differing Lorentz numbers. 

One first notices that Ohm's law @ leads to rewrite 
the equation ([5]) for the particle dynamics as 

X = -[u{X, t)-X] + h[X - u{X, t)]xB + rjj). (8) 

T t 

To get an impression of the relative importance of the 
three terms appearing on the right-hand side, Figure [5] 
depicts their root-mean-square values as a function of the 
Lorentz number. For small Lo the drag force dominates. 
Around Lo = 0.8 the electromagnetic part wins over the 
drag force. 



((U-V)/t) — 1— 
((u-V)xB/l) -X - 
(Tlj/I) 


„ -X 


.X 




* 

'■^ i 





0.5 1 1.5 2 2.5 3 

Lo 

FIG. 9: Rms-contributions of the different terms in ([8]l to the 
force of charged particles. 




-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 
(0 B 

FIG. 10: Probability density of the scalar product tu-B sam- 
pled along particle paths. Inset; skewness 5* of lj B. 



A quantitative measure of the preferential selection of 
certain sheets is provided by the skewness {{cu ■ B)^)/{{u! ■ 
^y^3/2 ^j^g distribution of the scalar product lu ■ B. 
As is shown in the inset of Fig. [10] the maximal skewness 
and hence the maximal preference for attractive sheets 
is attained for Lo ~ 1. There are two possible reasons 
for the decrease in the preferential selection of the sheets 
beyond this value. A first reason relies on a decreasing 
Larmor radius ( ~ 4?^ for Lo = 0.84), since the argument 
presented for the attraction and repulsion only holds for 
Larmor radii much larger than the Kolmogorov-scale. A 
second reason originates from a reduction of the degree 
of concentration at large Lorentz numbers, as will be dis- 
cussed in the following two subsections. 



B. Small-scale clustering 



A. Current sheet selection 

We first report some results on the effect of the charge, 
and thus of the Lorentz force acting on the particles, on 
their clustering properties. While neutral particles tend 
to cluster in the vicinity of any current /vortex sheet, 
charged particles are observed to concentrate preferen- 
tially on a subset among these sheets. Indeed, because 
of the Lorentz-force — u{X, t)] x B) exerted on the 
particles, there are attractive and repulsive sheets. In 
general the velocity of heavy particles is lower than the 
fluid velocity at the particle position. Hence, positively 
charged particles are attracted by sheets with u: ■ B > 0. 
Figure [TOl represents the probability density of this scalar 
product along the trajectories of neutral and charged par- 
ticles with the same Stokes number. One can observe 
that the distribution is tilted towards positive values. 



As for the neutral particles we measured the correla- 
tion dimension D2 (see Fig. fTTj) . In order to have a visible 
effect the Lorentz number has to exceed a certain thresh- 
old. Beyond this relatively flat region the charge leads to 
an increased concentration below the Kolmogorov scale, 
hence to a smaller D2. This is in relation to the prefer- 
ential selection of sheets discussed in the previous sub- 
section. The clustering is enhanced as the attraction ar- 
gument provides an additional mechanism for building 
concentrations. 

However, larger Lorentz numbers weaken the preferen- 
tial concentration significantly. This can be explained by 
the fact that the principal clustering process arises from 
the drag force. As the electromagnetic force wins over 
the drag force for large Lorentz numbers its influence is 
weakened and the correlation dimension exceeds that of 
the neutral particles, therefore reducing the small-scale 
clustering. 
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FIG. 11: Correlation dimension for ions 



C. Inertial-range clustering 

In the previous subsection we showed that small Lo 
increase the clustering at scales below the Kolmogorov 
scale. Turning to the implications of a finite Lorentz 
number for the inertial range of scales the qualitative 
findings are comparable. Fig. [12] confirms the tendency 
to an increased clustering for small charges also in the 
inertial range, with again a decrease at higher charges. 



1.05 




FIG. 12: Second order moments of the coarse grained density 
over a scale r = 18?7, as a function of the Lorentz number Lo, 
normalized by its value at io = 0. 

Note that the rescaling of the inertial-range distribu- 
tion as a function of the contraction factor F = {T]/r)St 
observed for neutral particles does not apply in the case 
of charged particles. The presence in that case of the 
Lorentz force, in addition to the viscous drag, does not 
allow an approach similar to that presented in Sect ion ITVl 



VI. CONCLUDING REMARKS 

We have presented results on the study of heavy parti- 
cle transport by MHD turbulent flows. Because of their 
inertia such particles were observed to display strong in- 
homogeneities in their spatial distribution. More pre- 
cisely, heavy particles tend to concentrate in the vicinity 
of the current sheets that are present in the conducting 
flow. Neutral (uncharged) particles approaches sheets, 
irrespectively of the nature of the latter, while charged 
particles, because of a combined effect of the Stokes drag 



and the Lorentz force, select those sheets where the scalar 
product between the magnetic field and the vorticity is 
of the same sign as their charge. 

The particle dynamics considered in this work is en- 
tirely passive, namely they have no feedback on the car- 
rier flow. To estimate their possible influence on the flow 
we computed the electric field Ep produced by the passive 
charged particles. This field is of course pronounced in 
regions with high concentrations. It is oriented in a way 
to act against the outer electric field Ej due to the charge 
in the fiow. To quantify these findings we measured the 
correlation C(r) = {Ep{x + r)Ef{x)) and observed that 
it was negative at all spatial scales. As seen from Fig. [T^ 
this correlation is observed to behave as a function of the 
scale like a power-law with exponent 1.7 in the dissipa- 
tive range and 0.6 in the inertial range. The origin of 
these power laws is still an open question which is clearly 
related to the concentration properties of the particles. 
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FIG. 13: Correlation of the flow electric field Ej and the elec- 
tric field induced by the particles E^: dashed (black online) 
line corresponds to a scaling index of 0.6, solid (red online) 
line to a scaling index of 1.7. 



It seems crucial to include electrical interactions be- 
tween particles as soon as the latters are not very diluted 
in the conducting fiow. Electric repulsive forces should 
clearly alter both the small-scale and the inertial-range 
properties of the particle spatial distribution but we still 
expect the preferential selection of current sheets to occur 
at inertial-range scales. 

Finally we would like to stress a possible application of 
the results presented here. Light particles, such as bubble 
in water, have now been used for a long time in order to 
give a direct visualization of the intense vortex filaments 
that are present in hydrodynamic turbulence. We have 
seen here that heavy particles in magneto-hydrodynamics 
tend to cluster in the vicinity of the current sheets that 
represent the most violent structures present in the tur- 
bulent motion of conducting flows. An idea could then 
be to use such particles in order to develop a new direct 
visualization technique of MHD flows that is expected to 
be very useful for grabbing qualitative and quantitative 
global information in turbulent dynamo experiments. 
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